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ABSTRACT 

The GJ876 system was among the earliest multi-planetary detections outside of the Solar System, 
and has long been known to harbor a resonant pair of giant planets. Subsequent characterization of 
the system revealed the presence of an additional Neptune mass object on an external orbit, locked 
in a three body Laplace mean motion resonance with the previously known planets. While this 
system is currently the only known extrasolar example of a Laplace resonance, it differs from the 
Galilean satellites in that the orbital motion of the planets is known to be chaotic. In this work, 
we present a simple perturbative model that illuminates the origins of stochasticity inherent to this 
system and derive analytic estimates of the Lyapunov time as well as the chaotic diffusion coefficient. 

We then address the formation of the multi-resonant structure within a protoplanetary disk and show 
that modest turbulent forcing in addition to dissipative effects is required to reproduce the observed 
chaotic configuration. Accordingly, this work places important constraints on the typical formation 
environments of planetary systems and informs the attributes of representative orbital architectures 
that arise from extended disk-driven evolution. 

1. INTRODUCTION 

The realization that planets form in gaseous protoplan¬ 
etary disks dates back to the collective works of Swe¬ 
denborg, Kant and Laplace. Despite several alterations, 
this nebular hypothesis has survived the test of time 
(Safronov 1969; Wetherill & Stewart 1989). Neverthe¬ 
less, a quantitative characterization of the consequences 
of planet-disk interactions (Kley & Nelson 2012) and the 
associated sculpting of orbital architectures of planetary 
systems (Morbidelli 2013) has only become an active field 
of research comparatively recently. 

Large-scale orbital migration of giant planets was first 
recognized as a theoretical possibility by Goldreich & 

Tremaine (1980), in an effort to quantify the orbital evo¬ 
lution of satellites embedded in circum-planetary disks. 

However, it wasn’t until the discovery of the first close-in 
extrasolar planets (Mayor & Queloz 1995) that this idea 
gained wide-ranging traction (Lin et al. 1996). Accord¬ 
ingly, over the last two decades, disk-driven migration 
has been repeatedly invoked as a unifying mechanism to 
explain the various orbital properties of extrasolar giant 
planets (Lee & Peale 2002; Grida et al. 2007; Batygin 
2012) as well as low-mass compact multi-planetary sys¬ 
tems alike (Papaloizou et al. 2007). 

Qualitatively distinct modes of migration^, character¬ 
ized by different timescales and directions, can arise de¬ 
pending on the physical properties and structure of the 
disk as well as the embedded planet (Ward 1997; Grida et 
al. 2007; Paardekooper & Papaloizou 2008, 2009; Bitsch 
& Kley 2011; Bitsch et al. 2013). Because of this intrinsic 
diversity in physical behavior, simultaneous migration of 
multiple planets residing within the same disk may cause 

^ Under simplifying assumptions, the two most-commonly 
quoted categories of migration are that associated with the vis¬ 
cous evolution of the disk (characteristic of giant planets that are 
able to clear out substantial gaps in their orbital neighborhood) 
and that facilitated by resonant interactions (characteristic of low- 
mass planets that remain immersed in the gas). 


the orbits to approach each other. 

Convergent orbital evolution can result in capture of 
objects into mean-motion resonances (orbital states char¬ 
acterized by rational period ratios) and evidence for this 
process is plentiful throughout the satellite population 
of the Solar System (Goldreich 1965; Yoder 1973, 1979; 
Greenberg 1977; Peale 1976, 1986; Henrard 1982; Hen- 
rard & Lemaitre 1983; Malhotra 1990). Moreover, the 
existence of a substantial number of (near-)resonant ex¬ 
oplanet ary systems (see Wright et al. 2011; Batygin & 
Morbidelli 2013) suggest that resonant locking is not lim¬ 
ited to satellites and is also active among planets (Masset 
& Snellgrove 2001; Morbidelli et al. 2007). 

Despite their non-negligible count, resonant systems 
comprise a minority within the current observational ag¬ 
gregate (Wright et al. 2011; see also Fabrycky et al. 
2014). Instead, giant planets often reside on dynami¬ 
cally excited orbits that are believed to be a result of 
planet-planet scattering (Rasio & Ford 1996; Juric & 
Tremaine 2008; Raymond et al. 2009). The connection 
between processes that occur within the nebular epochs 
of planetary systems and the onset of large-scale or¬ 
bital instabilities responsible for sculpting the observed 
semi-major axis - eccentricity distribution are at present 
poorly understood (Lega et al. 2013). Nevertheless, it 
is entirely plausible that planet-planet scattering origi¬ 
nates within compact systems assembled by disk-driven 
migration (Beauge & Nesvorny 2012; Morbidelli 2013). 

It is worthwhile to note that the presently favored view 
of the Solar System’s early dynamical history (see Mor¬ 
bidelli et al. 2008 for a review) is consistent with the pic¬ 
ture delineated above. That is, an evolutionary sequence 
where the giant planets emerged from the Solar nebula 
in a multi-resonant configuration that subsequently be¬ 
came unstable, causing the planets to scatter onto their 
current orbits is broadly consistent with the available 
observations (Tsiganis et al. 2005; Morbidelli et al. 2007; 
Levison et al. 2008; Batygin et al. 2011; Nesvorny & Mor- 
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bidelli 2012). 

Numerous insights into the physical processes that op¬ 
erate within birth environments of planetary systems 
could be obtained if the orbital states at the time of 
nebular dispersion could be inferred. In practice, this is 
difficult to do for dynamically relaxed systems because 
planet-planet scattering is a highly stochastic “forward” 
process that tends to erase all memory of the system’s ini¬ 
tial state. In other words, the characterization of the vast 
majority of the observational sample can yield only lim¬ 
ited information about the primordial nature of the asso¬ 
ciated planetary systems because their dynamical histo¬ 
ries have been chaotically eliminated. Even in the Solar 
System, where the modeling efforts have enjoyed a mul¬ 
titude of observational restrictions, strong degeneracies 
with respect to the initial state persist (Batygin & Brown 
2010; Nesvorny 2011; Batygin et al. 2012). 

On the contrary, resonant planetary systems (which 
have managed to escape the onset of instability^) may in 
fact yield tangible constraints on the environment within 
which they formed. As such, they constitute exception¬ 
ally high value targets for theoretical inquiries. 

Arguably the most exotic resonant exoplanetary sys¬ 
tem detected to date is GJ876 (Marcy et al. 2001; Rivera 
et al. 2005), and this remarkable collection of planets 
will be the primary subject of this paper’s study. With 
three objects locked in resonance, GJ876 comprises the 
only well-characterized extrasolar multi-resonant chain, 
although additional multi-resonant systems exist within 
the Kepler data set (e.g. Kepler-16, Kepler-79, KOI-730; 
Fabrycky et al. 2014). While this system exhibits some 
similarity to the Galilean satellites (Peale 1976, 1986), 
it differs from the lo-Europa-Ganemede system in a cru¬ 
cial aspect: the resonant arguments of GJ876 exhibit 
vigorous, yet bounded (on multi-Gyr timescales) chaos 
(Rivera et al. 2010; Marti et al. 2013). 

The architecture of GJ876 is fully consistent with the 
picture of conventional disk-driven migration (Lee & 
Peale 2002; Crida et al. 2008). As a result, the dynami¬ 
cal characterization of GJ876 provides a rare window into 
the description of orbital properties of planetary systems, 
as they appear when they emerge from their natal disks. 
With this notion in mind, performing a theoretical anal¬ 
ysis of the system’s dynamical state with an eye towards 
obtaining some insights into the nature of the disk from 
which this system formed is the goal of this work. 

The paper is structured as follows. In section 2, we 
examine the qualitative features of the system’s orbital 
architecture and setup the basis for theoretical analysis. 
In section 3, we construct a perturbative model for the 
resonant dynamics of the system and elucidate the ori¬ 
gins of chaotic motion in an analytical fashion. In section 
4, we consider the assembly of the multi-resonant chain 
in a dissipative and turbulent protoplanetary nebula. We 
conclude and discuss our results in section 5. 

2. THE PHYSICAL SETUP OF THE CALCULATION 

The observational saga of the GJ876 system effectively 
spans the entire active history of exoplanetary science. 
The initial detection of a ^ 2Mjup giant planet “b” in a 
61 day orbit dates back to the infancy of large-scale dedi¬ 
cated radial velocity surveys (Delfosse et al. 1998; Marcy 

^ See Raymond et al. (2008) for an alternative view-point. 


TABLE 1 

Adopted Orbital Fit of the G,I876 System 



M (Mq) 

a (AU) 

e 

M (deg) 

ro (deg) 

★ 

0.334 

- 

- 

- 

- 

d 

2.051 X 10“® 

0.0208 

0.207 

355 

234 

c 

6.820 X 10-* 

0.1296 

0.256 

294.59 

48.76 

b 

2.173 X 10“® 

0.2083 

0.032 

325.7 

50.3 

e 

4.385 X 10“® 

0.3343 

0.055 

335 

239 


et al. 1998). The discovery of a ~ 0.7Mjup companion 
“c” on a 30 day orbit followed shortly thereafter (Marcy 
et al. 2001), rendering the GJ876 “c-b” pair the first 
mean motion resonance to be identified outside the Solar 
System. Taking advantage of the observational imprint 
of resonant coupling, Laughlin & Ghambers (2001) and 
Rivera & Lissauer (2001) presented independent analy¬ 
ses of the radial velocity data that accounted for planet- 
planet interactions and were able to break the sin(i) de¬ 
generacy inherent to radial velocity detections, deriving 
the system’s inclination with respect to the line of sight 
of J ~ 50 deg. Concurrently, the signal of planet “b” was 
confirmed astrometrically by Benedict et al. (2002) (see 
also Bean & Seifahrt 2009). 

Detections within the system continued, as an addi¬ 
tional close-in ~ 7.5M0 planet “d”, residing on a 2 day 
orbit was announced by Rivera et al. (2005). Under the 
assumption of a 3-planet system, Correia et al. (2010) 
re-analyzed the available data and with extensive mod¬ 
eling derived a mutual inclination between the resonant 
planets “b” and “c” of Ai ~ I deg. Moreover, the study 
of Correia et al. (2010) confirmed the existence of planet 
“d” and strongly hinted at the eccentric nature of its 
orbit (see also Baluev 2011). 

The latest advancement in the observational charac¬ 
terization of GJ876 arose from the work of Rivera et al. 
(2010), who uncovered yet another resonant ~ IbM^ 
planet “e,” occupying a 124 day orbit. Further dynam¬ 
ical analysis revised the system inclination to i ~ 60 
deg and more importantly, showed that the evolution of 
the multi-resonant configuration undergoes bounded, yet 
chaotic variations. 

Over the last decade and a half, the orbital state of 
the system and its origin have been studied by a sub¬ 
stantial number of authors, employing a wide variety of 
methods. Specifically, in addition to the studies quoted 
above, orbital characterization and stability have been 
explored by: Jones et al. (2001); Kinoshita & Nakai 
(2001); Gozdziewski et al. (2002); Ji et al. (2002); Zhou 
& Sun (2003); Haghighipour et al. (2003); Beauge & 
Michtchenko (2003); Veras (2007); Baluev (2011) and 
Marti et al. (2013). Meanwhile the assembly of the par¬ 
ticular resonant configuration has been simulated and 
studied by: Snellgrove et al. (2001); Murray et al. (2002); 
Lee & Peale (2002); Thommes & Lissauer (2003); Beauge 
et al. (2003, 2006); Kley et al. (2004, 2005); Lee (2004); 
Adams et al. (2008); Crida et al. (2008) and Lee & 
Thommes (2009). Additionally, Gerlach & Haghighipour 
(2012) examined the possibility of extra bodies being 
locked in the multi-resonant configuration. 

Thanks to the long time-span covered by the radial 
velocity observations and the aforementioned fitting ef¬ 
forts, the orbital properties of planets “b” and “c” are 
well constrained (Correia et al. 2010). However, sub- 
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Fig. 1.— Dynamical evolution of the isolated “c-b” resonant pair. The left panels (A and D) show semi-major axis evolution, while the 
right panels (B,C,E and F) show eccentricity and critical argument evolution, expressed in terms of scaled canonical cartesian variables. 
The top and bottom rows corresponds to planets “c” and “b” respectively. In each plot, blue curves show osculating elements obtained from 
direct A^-body integration. The red curves denote numerically averaged output, where the high-frequency component has been removed 
by Fourier analysis (Laskar 1993a; Morbidelli 1993). Treating the outermost planet as a massless test-particle renders the evolution of the 
massive resonant pair quasi-periodic. 


stantial uncertainties exist in the orbital solutions of the 
lower mass components “d” and “e” (Rivera et al. 2010). 

The desired precision of the knowledge of the dynam¬ 
ical state of a system is largely dictated by the purpose 
of the calculation one wishes to perform. Here, our aim 
is not to create an ephemeris for GJ876, but rather to 
shed light on the origins of chaotic motion within the 
multi-resonant system and place rough constraints on 
the properties of the nebula within which the system was 
born. To this end, we note that once a 4-planet system 
is adopted, chaos (highlighted by the stochastic evolu¬ 
tion of planet “e”) is more or less ubiquitous throughout 
the parameter space restricted by the data (Marti et al. 
2013). Consequently, for the purposes of this work, we 
shall simply adopt the best-fit co-planar orbital solution 
of (Rivera et al. 2010) at face-value, keeping in mind 
that quantitatively different evolutions stemming from 
nearby initial conditions can be equally representative of 
the system’s dynamical behavior. To this end, we fur¬ 
ther remark that while the orbital fit of Rivera et al. 
(2010) clearly favors a chaotic solution, there may exist 
quasi-periodic islands in phase space that reside within 
parameter space covered by observational uncertainties. 
For completeness, the adopted orbital solution is pre¬ 
sented in Table (1). Note that in this work, we have 
adopted a host stellar mass estimate of Correia et al. 
(2010), which differs from that of Rivera et al. (2010) by 
^ 4%, meaning that our results will differ from previous 
works on a detailed level. 

For the purposes of understanding the evolution of the 
“c-b-e” multi-resonant chain in the simplest terms pos¬ 
sible, the perturbations arising from the close-in planet 
“d” can be neglected^. Broadly speaking, this is justified 

^ Admittedly, doing so removes some of the high-order features 


by the fact that while each member of the chain inter¬ 
acts with its nearest neighbor(s) via first order resonant 
terms that scale as oc e (where e is the eccentricity), 
the averaged gravitational coupling between planet “d” 
and the multi-resonant system is secular in nature and to 
leading order scales as oc (Murray & Dermott 1999). 
To confirm this argument, we numerically integrated the 
V-body evolution of the system using the mercury6 soft¬ 
ware package (Chambers 1999) with and without planet 
“d” and found no meaningful differences in the obtained 
solutions. 

It is further noteworthy that the mass of the outer¬ 
most planet is more than an order of magnitude smaller 
than the two inner members of the chain. It is thus 
tempting to further simplify the system and treat planet 
“e” as a test particle, subject to perturbations arising 
from the “c-b” resonant pair. Numerical examination of 
such a setup reveals that qualitative aspects of the dy¬ 
namical evolution of planet “e” are largely insensitive to 
this assumption. On the contrary, under this simplifi¬ 
cation the evolution of the massive resonant pair ceases 
to be chaotic and instead exhibits quasi-periodic motion 
(Beauge et al. 2003; Correia et al. 2010). 

Figure (1) shows the results of the simplified integra¬ 
tion where planet “e” is treated as a test particle, and 
planet “d” is neglected. Evidently, the stochastic fea¬ 
tures of the multi-resonant chain arise entirely as a re¬ 
sult of the perturbations of planets “c” and “b” onto 
planet “e,” and the chaotic diffusion of the massive pair 
is communicated exclusively via a back-reaction. This at¬ 
tribute has important consequences for the purposes of 
constructing a simple model aimed at elucidating the dy- 

of the system’s dynamical portrait. However, detailed scrutiny is 
unimportant to the problem at hand. 


































4 


- TOb - Azu (degrees) = ebcos(rob) -N-body - analytical 



kb = eb sin(tnb) 



Fig. 2. — Dynamical evolution of planet “b” in the periodic approximation. The time evolution of the longitude of perihelion of planet 
“b” is shown with blue points in panel A. Concurrently, the apsidal difference between planets “b” and “c” is shown as a red curve in panel 
A. Panels B and C depict the evolution of the scaled canonical cartesian coordinates corresponding to planet “b.” As in Figure (1), blue 
curves represent the time-series obtained from direct numerical integration. Meanwhile, the over plotted red curves show the analytical 
parameterization given by equation (2). 


namical structure that underlies chaotic motion. Specif¬ 
ically, in this work we will take advantage of this charac¬ 
teristic and derive the chaotic properties of the system by 
considering the evolution of a massless planet “e,” sub¬ 
ject to periodic perturbations from planet “b,” whose 
orbit is in turn modulated by interactions with planet 
“c.” 

Panel A of Figure (2) shows the numerically obtained 
time evolution of the longitude of perihelion of planet 
“b” (blue) as well as difference between the longitudes of 
perihelia of planets “c” and “b” (red). Clearly, the “c-b” 
resonant pair resides near a dynamical equilibrium char¬ 
acterized by alignment and co-precession of the apsidal 
lines of the orbits (Laughlin & Chambers 2001; Lee & 
Peale 2002). Such a state is actually quite peculiar for 
resonant orbits and is only possible when the orbital ec¬ 
centricities are not small. Indeed, the classical first-order 
expansion of the resonant Hamiltonian (Leverrier 1855; 
Ellis & Murray 2000) does not capture the existence of 
this fixed point. To this end, Beauge et al. (2003) de¬ 
veloped an alternative expansion of the planetary three- 
body Hamiltonian and showed that the resulting pertur¬ 
bative (first-principles) solution matches the numerically 
computed evolution well. Taking a somewhat alternative 
approach, Correia et al. (2010) utilized synthetic pertur¬ 
bation theory to construct a Lagrange-Laplace like peri¬ 
odic solution for the resonant pair which also shows ex¬ 
cellent agreement with iV-body calculations. Because the 
dynamics of the “c-b” resonant pair have been studied 
extensively, here we shall not duplicate the published re¬ 
sults and instead refer the interested reader to the afore¬ 
mentioned studies. 

As already briefly mentioned above, in this work we 
shall examine the consequences of gravitational excita¬ 
tion of planet “e” facilitated by planet “b.” In order to 
perform this analysis within a perturbative framework, 
we must first delineate an approximate functional form 
for the dynamical behavior of planet “b.” Let us begin 
by defining the following scaled cartesian canonical coor¬ 
dinates: 

h = ecos(tJ7) 

k = esm{zu), ( 1 ) 

where w is the longitude of perihelion. Following Cor¬ 
reia et al. ( 2010 ), we note that in the {h,k) plane the 
trajectory of planet “b” is well described by a circle of 
radius 6 that is off-set from the origin by e (see also the 


discussion on free and forced elements in Ch. 7 of Mur¬ 
ray & Dermott 1999). Accordingly, we parameterize the 
evolution of planet “b” in the following manner: 

hh = e + S cos{gt + vuq) 

kh = S sm{gt + vjq). (2) 

In the above expression, g = —0.6706 rad/yr is the 
(retrograde) apsidal precession rate and wq is the phase 
offset, while the constants are set to e = 0.004 and 
6 = 0.035. Panels B and C in Figure (2) depict a compar¬ 
ison between results obtained with an IV-body simulation 
and the analytical prescription (2). Needless to say, the 
observed agreement is excellent. 

It is further noteworthy that the numerically averaged 
semi-major axis of planet “b” does not deviate from its 
nominal value, [a]b by more than a few parts in a thou¬ 
sand. Thus, for the purposes of the following calculation 
we may readily neglect the semi-major axis evolution 
all together and set Ob = [a]b. With a simple analyti¬ 
cal model for the dynamical evolution of the perturbing 
planet at hand, we are now in a position to perform a 
perturbative analysis of the system’s chaotic behavior. 

3. ORIGINS OF CHAOTIC MOTION 

Before we proceed, let us begin by defining restricted 
Poincare action-angle variables in terms of standard or¬ 
bital elements: 

A = X = Ai + uj 

F = A(1 - a/I - e2) ~ [A]eV2 7 = -w (3) 

where Q is the gravitational constant, M* is the stellar 
mass, M is the mean anomaly and the nominal action 
[A] is evaluated at [a] = 2^/^[a]b. For the remainder 
of the paper, quantities not labeled with an index shall 
correspond to planet “e.” 

3.1. Resonant Perturbation Theory 

Generally, unlike secular perturbations (see e.g. Laskar 
1996), mean-motion resonances modulate both the ec¬ 
centricity and semi-major axis. Hence, Keplerian motion 
cannot be averaged over, and the corresponding part of 
the Hamiltonian must be retained: 

^ 

«-kep — „ . 2 


( 4 ) 
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Fig. 3.— Chaotic evolution of the outermost planet, “e.” Panel A shows the eccentricity as a function of time, while panels B and C show 
the evolution of the critical arguments corresponding to the first-order (2:1) mean motion resonance. Across all panels, the blue curves 
represent time series obtained by direct A^-body integration, the black curves denote solutions arising from the autonomous two degree of 
freedom perturbative model (Hamiltonian 22), and the green curves show the evolution described by the simplified non-autonomous one 
degree of freedom model (Hamiltonian 28). Note that all solutions track each-other well initially but lose coherence after a few decades of 
evolution. 


Expanding "Hkep around the nominal 2:1 resonant semi¬ 
major axis quoted above, equation (4) takes the form: 


= -SS + ^(A- [A]) - [A])^ 


2[A]2 [A]3 


[A]2 


4[n]A--[h]A , 


( 5 ) 


where [n] = ^5M*/[a]3 is the mean motion and [h] = 
[n]/[A]. The first term on the second line of equation (5) 
is constant and can thus be dropped from the Hamilto¬ 
nian. 

Upon averaging out short-periodic terms (see Ch. 2 
of Morbidelli 2002), the component of the Hamiltonian 
that governs first order (in e) resonant planet-planet in¬ 
teractions (within the framework of the elliptic restricted 
three body problem) reads (Murray & Dermott i999): 


'Hres = - (^/re^-’eb COs(2A - Afa - Wb) 

/ 2r \ 

+/reiW pq cos(2A - Ab-k 7)j , (6) 

where fres = —1.1905 and fres = 0.4284 are interaction 
coefficients that depend exclusively on the semi-major 
axis ratio (Murray & Dermott 1999). The presence of 
Cb and vjh in the Hamiltonian (6) gives rise to explicit 
time-dependence in expression (6). Thus, the quoted 
model constitutes a non-autonomous dynamical system 
with two degrees of freedom. 

The Hamiltonian can be made autonomous by extend¬ 
ing the phase-space to three degrees of freedom (Mor¬ 
bidelli 2002): 

■H = Hkep +'Hres + T, (7) 

where T is the newly introduced action conjugate to time, 
t. Accordingly, substituting the solution (2) into equa¬ 
tion (6), we obtain the following expression: 

n = A[n]A-^[h]A^+T-^ 

X ((e + ^ cos{gt + (po)) cos(2A - n^t - Abo) 

-I- Ssin{gt — wq) sin(2A — ribt — Abo)) 




cos(2A - riht - Abo + 7 ) 


( 8 ) 


where Abo is an initial phase of planet “b,” and Uh is its 
mean motion. 

Equation (8) is rather cumbersome and can be made 
more succinct. First, let us define the constants: 


a = — 



GMb fiel 

H yiAj' 


(9) 


Next, let us perform a canonical transformation of co¬ 
ordinates, arising from the following type-2 generating 
function (Goldstein 1950): 

T 2 = (2A — ribt — Abo)0 + ( 7 )*^* + (10) 


An application of the transformation equations yields the 


new action-angle variables: 
0 = A/2 

$ = r 

S = T -I- rib© 


0 = 2A — Tibt — Abo 
(/) = 7 

e = (11) 


In terms of the new coordinates, the Hamiltonian reads: 

% = 8[n]0 — 6[/i]0^ -I- ae cos{6) -\- a5 cos(0) 

X cos{g^ -I- ifo) + aS sin(0) sin(55 -|- ipo) 

+ cos{9 + (jj) — nb& + ^. (12) 


Although somewhat less unwieldily than equation (8), 
the Hamiltonian (12) is still characterized by three de¬ 
grees of freedom, precluding straightforward analytical 
treatment. To remedy this issue, let us define canoni¬ 
cal cartesian coordinates related to the ($, </) degree of 
freedom: 

x = V^cos{(j)) y = •\/^sin((/)). (13) 

After some algebraic manipulation, the Hamiltonian ob¬ 
tains the form: 


H = (8[n] — nb)0 — 6[/i]0^ -|- aS cos{9 — g^ — po) 

+ ae cos{9) + j3x cos{9) — /3y sin(d) -|- S. (14) 
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To reduce the number of harmonics present in "H, we 
follow Henrard et al. (1986); Wisdom (1986) and define 
the following canonical translation"*: 

x = x+-e y = y- (15) 

This transformation morphs the first and second terms 
on the second line of equation (14) into a single term. Ac¬ 
cordingly, upon defining implicit action-angle variables 

X = cos{ip) y = sin^-ip), (16) 

the Hamiltonian takes on a form characterized by only 
two harmonics: 

% = (8[n] — nb)0 — 6[/i]0^ -|- aS cos(0 — — ^Po) 

+ cos{9 + tjj) + E. (17) 

From here, the reduction of the Hamiltonian to a two- 
degree of freedom system is straightforward. Particularly 
as in equation (10), define a contact transformation aris¬ 
ing from a type-2 generating function: 

T2 = {0-9^-Vo)W+{9 + ij)Z + iOIC. (18) 

The transformation equations yield: 

W — Q — Z w = 9 — — ipo 

Z = 'I' z = 9 + tjj 

IC = E + gW K = C (19) 


Accordingly, the Hamiltonian is rewritten as follows: 

— gW + a6 cos{w) + l3u. (22) 

The corresponding equations of motion take the form: 

—— = ao sm(r(;) 
at 

J=^ + 6[n]!i-12[%(yV+A^) (23) 


After numerous variable changes, it is useful to explic¬ 
itly relate the final variables (19,21) to Keplerian orbital 
elements. Working back through the transformations de¬ 
lineated above, it can be shown that the orbital semi¬ 
major axis and eccentricity are expressed as follows: 


a = 


e = 


{2W + {u-e{alP)) +v^) 


2)2 


GM^ 


{u — e{a/P)y + 




(24) 


Noting that the resonant condition implies that rib = 2[n] 
the Hamiltonian takes the form: 

U = 6 [n](>V + Z)- 6[h]{W + Zf - gW 

+ aS cos{w) + pV^ cos{z) + K.. (20) 

Because the angle k (which corresponds to time) is now 
absent from the Hamiltonian, the action /C is a constant 
of motion and can thus be dropped from equation ( 20 ). 

The reduction of the Hamiltonian to an autonomous 
two-degree of freedom system is now complete. Quali- 
tatiely, the Hamiltonian (20) represents a system of two 
momentum-coupled pendulums. It is noteworthy that 
in isolation, the {Z,z) degree of freedom possesses the 
D’Alembert characteristic and thus has the dynamics 
corresponding to the second fundamental model for res¬ 
onance (Henrard & Lemaitre 1983) whereas the isolated 
(yy, w) degree of freedom has the phase-space structure 
of a simple pendulum (Chirikov 1979). 

The extent to which the perturbative model governed 
by Hamiltonian (20) successfully captures the dynamical 
behavior of GJ876 “e” can be tested numerically. How¬ 
ever prior to applying Hamilton’s equations, note that 
in terms of action-angle coordinates, Hamiltonian (20) 
possesses a coordinate singularity at 2^ = 0. This nui¬ 
sance can be circumvented by transforming to canonical 
cartesian coordinates: 

u = V^cos{z) v = V^sm{z). ( 21 ) 

Note that for the non-restricted (i.e. planetary) three-body 
problem, there exists a corresponding canonical rotation that re¬ 
duces the first-order resonant Hamiltonian to an integrable one 
(Sessin & Ferraz-Mello 1984; Wisdom 1986; Batygin & Morbidelli 
2013; Deck et al. 2013; Delisle et al. 2014). 


Meanwhile, the angles present in the original formulation 
of the Hamiltonian ( 6 ) are related to the variables (19) 
in a unembellished way: 

w = (2A — Ab — vj\,) 

z = arctan ~ (2A — Ab -I- 7 ). (25) 

The latter equality is inexact, but is nevertheless a good 
approximation because e <C (5. In the same spirit, the 
following relationships approximately hold: 

u ~ ecos(z) 


u ~ e sin(z) 



Note that the last of the above expressions is closely re¬ 
lated to the well-known Tisserand parameter. 

Orbital evolution obtained by numerical integration of 
equations (23) is shown with black curves in Figure (3), 
where the results of an fV-body simulation are also de¬ 
picted with blue lines. Evidently, the system described 
by Hamiltonian (22) provides an excellent perturbative 
representation of the real system, meaning that the dy¬ 
namical behavior within the chaotic layer is well captured 
by a first-order expansion of the disturbing function. It 
is noteworthy that the semi-analytical solution initially 
tracks the A^-body solution, but the two evolutionary se¬ 
quences lose coherence after ~ 50 years. This is indica¬ 
tive of a Lyapunov time that is a factor of a few shorter 
than 50 years®. This places GJ876 into the same cate- 

® This is in some contrast with a 10^ — 10^ year Lyapunov time 
reported by Rivera et al. (2010). 
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Poincare surface of section 
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Fig. 4.— Surfaces of section corresponding to the autonomous two degree of freedom Hamiltonian (22). Panel A shows level curves of 
the Hamiltonian on a double section, and the rest of the panels depict Poincare surfaces of section with respect to the (W, w) degree of 
freedom, at various levels of "H. The nominal resonant semi-major axis is shown with a red curve in panel A. In the other panels, the 
variables are scaled such that the radial distance away from the origin approximately corresponds to the eccentricity. On the Poincare 
surfaces of section, the chaotic sea is shown with black points, while quasi-periodic trajectories are shown with purple, brown, and orange 
points. Note that chaotic excisions are limited by the conservation of "H, which give rise to inadmissible regions in phase-space. 


gory of rapidly chaotic systems as Kepler-56 (Deck et al. 

2012 ). 

Despite rapid chaotic diffusion, V-body calculations 
reported by Rivera et al. (2010) suggest that the system 
is stable on multi-Gyr timescales. In other words, the 
chaos exhibited by the multi-resonant system is simulta¬ 
neously vigorous and bounded. Both of these character¬ 
istics can be qualitatively understood within the context 
of Hamiltonian (22) by examination of surfaces of sec¬ 
tion. 

As a starting point, let us examine the evolution of 
the angles depicted in Figure (3). Note that the angle 
(2Ab — A — zuh) librates around 0 while (2Ab — A -I- 7 ) 
switches between libration around tt and circulation. 
Preliminary intuition about the amplitude of orbital ex¬ 
cursions can be obtained by constructing a double sur¬ 
face of section of the Hamiltonian. In accord with the 
evolution of angles denoted in Figure (3), we fix w = 0 
and u = 0 and plot level curves of T-L on panel A of 
Figure (4). Suitably, on the x-axis, we plot the scaled 
action 2W/[A] while on the y-axis we plot the scaled 
momentum u/-^/[Aj. Because the Hamiltonian (22) is 
autonomous, any (potentially chaotic) trajectory it de¬ 
scribes will be constrained to map onto a corresponding 
level curve (given by the value of TL) every time the or¬ 
bital state crosses the section condition. 

A red curve corresponding to the nominal semi-major 
axis a = [a] is also shown on this panel. Crudely speak¬ 
ing, the horizontal deviation away from the nominal 
resonance curve on the double section is indicative of 


changes in semi-major axis whereas vertical deviation 
corresponds to eccentricity modulation. Thus, the area 
occupied by a given constant H curve on the double sec¬ 
tion serves as a proxy for the amplitude of orbital vari¬ 
ations associated with resonant motion. This point is 
of some importance to understanding the long-term sta¬ 
bility of the system. Particularly, this simple analysis 
suggests that to the extent that Hamiltonian (22) pro¬ 
vides an adequate representation of the dynamics, the 
conservation of K restricts the maximal deviation from 
nominal resonance of the trajectory, no matter how vig¬ 
orously chaotic it may be. Qualitatively, this explains 
how a rapidly stochastic system such as GJ876 can re¬ 
main stable on multi-Gyr timescales. 

In addition to delineating energy contours on a double 
section, it is further useful to explore the actual dynami¬ 
cal behavior of the system on the contours (i.e. provided 
specific values of V.). To this end, we have constructed 
Poincare surfaces of section for the five highlighted curves 
shown on the double section, labeled Hi...'H 5 . The sec¬ 
tions are taken with respect to ic = 0 at dw/dt < 0 
and depict the {u,v) degree of freedom. The surfaces 
are shown as panels B-F on Figure (4) and are labeled 
according to the level of H they represent. The value of 
"H that corresponds to the initial condition depicted in 
Table (1) is and is shown on panel D. 

Quasi-periodic trajectories, shown as purple and brown 
curves dominate the surfaces of section that correspond 
to energy levels that yield the most limited orbital vari¬ 
ations (i.e. Tii and H 5 , shown on panels E and F re- 
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stroboscopic surface of section 
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Fig. 5.— Properties of the simplified model governed by Hamiltonian (28). Panel A depicts the evolution of the action W, where the 
color scheme is the same as that employed in Figure (3). Note that the evolution corresponding to the non-autonomous system is simply 
that dictated by equation (27). Panel B shows a stroboscopic surface of section with the same section conditions as that employed in Figure 
(4). The sweeping of the separatrix is shown by depiction of its sequential shape at different times. The maximal extent of the hyperbolic 
equilibrium point (tihyp) is also explicitly labeled. In order to highlight the non-adiabatic nature of the system, panel C shows the adiabatic 
limit (where the modulation frequency cr is reduced by a factor of 10) of the equation of motion arising from Hamiltonian (28). Note that 
reducing a significantly introduces quasi-periodic resonant trajectories. 


spectively). On the contrary, the remaining surfaces of 
section feature substantial chaotic seas, shown with black 
points. The location of a given point on the surface of 
section is related to the orbital eccentricity (evaluated 
at the section conditions) through expression (24). Ac¬ 
cordingly, the range of chaotic eccentricity excursions at 
a given energy level can be gathered from examining the 
extent of the phase-space occupied by the irregular re¬ 
gion. Note further that in all cases, the dynamics resides 
on a well-defined domain, which is restricted by the con¬ 
servation of H and the requirement for the actions to 
have a null imaginary component. 

The physical characters of the periodic points that re¬ 
side at the centers of the nested quasi-periodic trajecto¬ 
ries is not uniform across the plotted energy levels. Par¬ 
ticularly, on panels B and C, these orbits correspond to 
limit cycles that are characterized by rapid circulation of 
w and a modulated circulation of z, such that the ma¬ 
jority of time spent on the circulation cycle remains in 
the vicinity of z = tt. Conversely, fixed points in panels 
D, E, and F are true equilibria that are characterized by 
stationary evolution of the angles at w = 0 and z = tt. 

Despite these differences, as long as the system is initi¬ 
ated within the extensive chaotic region (irrespective of 
whether the value of the Hamiltonian is set to Hi, H 2 , or 
H 3 ), the qualitative behavior of the dynamics resembles 
that observed in Figure (3). In light of this, it is tempt¬ 
ing to obtain an estimate for the characteristic rate of 
chaotic decoherence of the system, that is independent 
of the exact value of H. The simplest approach to such 
an estimate requires further reduction of the model. 

3.2. A Simplified Model 

Recall that the only assumptions inherent to the 
Hamiltonian (22) are the truncation of the disturbing 
function at first order in e (Murray & Dermott 1999) and 
the adoption of equations (2) for the dynamical evolution 
of the perturbing planet. The subsequent conversion of 
the Hamiltonian to an autonomous two degree of freedom 
system was possible because of the reducing transforma¬ 
tion (15) (Henrard et al. 1986; Wisdom 1986). Because 
the Hamiltonian (22) cannot be simplified further with 
canonical transformations, in order to convert the sys¬ 
tem into a more tractable non-autonomous one degree of 


freedom system (see e.g. Timofeev 1978; Escande 1985), 
we must prescribe a functional form to one of the degrees 
of freedom. 

An examination of panels B and C of Figure (3) shows 
that while the oscillations of the angle w are approxi¬ 
mately limited to the range —7r/2 ^ w < Tr/ 2 , the angle z 
recurrently transitions between libration and circulation. 
This implies that the repeated separatrix crossing associ¬ 
ated with the (u, v) degree of freedom acts as the primary 
driver of stochasticity (Lichtenberg & Lieberman 1983). 
Moreover, recall that the equations of motion (23) dictate 
that the interactions between the degrees of freedom are 
facilitated exclusively by action coupling. Consequently, 
for further reduction of the model, it is sensible to pa¬ 
rameterize the time-evolution of W. 

In absence of coupling, Hamiltonian (22) describes sim¬ 
ple pendulum-like dynamics for the (W, w) degree of free¬ 
dom. Under the assumption of small-amplitude libration 
of w (which allows one to approximate the dynamics of 
a pendulum with that of a harmonic oscillator), W will 
exhibit sinusoidal variations (Goldstein 1950; Morbidelli 
2002). Accordingly, let us adopt the following functional 
form: 

yV=^^{r] + fiCOs{at)). (27) 

Empirically, the newly introduced constants and fre¬ 
quency are set to p = 0.99, /i = 0.025, and a = 
27r/(14.72) rad/year, to provide a suitable match to the 
numerical calculations. To this end, panel A of Figure (5) 
shows a comparison between the evolution of W, com¬ 
puted within the framework of an A^-body simulation 
(blue), numerical integration of the perturbative model 
governed by Hamiltonian (22) (black) and the prescrip¬ 
tion (27) (green). For consistency, a sinusoidal wave with 
the same frequency (corresponding to the envisioned evo¬ 
lution of w, given the approximation 27) is also depicted 
on panel B of Figure (3). We note that even though the 
adopted parameterization does not respect the physical 
requirement for the quantity 2yV’/A to not exceed unity 
(see equation 26), it suffices for the purposes of an illus¬ 
trative model we aim to construct. 

Adopting equation (27), equations of motion (23) can 
again be integrated numerically to yield an approximate 
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dynamical evolution for the (u, v) degree of freedom. The 
resulting solutions for the eccentricity as well as the criti¬ 
cal angle (2A —Ab+y) are over-plotted as green curves in 
panels A and C of Figure (3) respectively. Although the 
shown solutions lose coherence (as chaotic systems must) 
after a few tens of years of evolution, the simplified model 
clearly captures the stochastic behavior exhibited by the 
system in a satisfactory manner. Thus, it can be sensi¬ 
bly employed to further characterize the dynamics in a 
rudimentary fashion. 

To understand the origins of chaos observed in the 
numerical solutions, note that under the assumption 
(27), the equations of motion (23) correspond to a non- 
autonomous one degree of freedom Hamiltonian: 

* = 6[nl(F^)+;3„ 

- (^^ + ^ 

This Hamiltonian possesses a single critical curve that 
sweeps across a region of phase-space with every circu¬ 
lation of the angle at (Henrard 1982; Cary et al. 1986). 
Accordingly, repeated encounters with the critical curve 
(which comprises an infinite-period orbit) render the mo¬ 
tion on the separatrix-swept region of phase-space irreg¬ 
ular (Chirikov 1979; Wisdom 1985). 

Panel B of Figure (5) depicts the critical curve of 
Hamiltonian (28) at various phases of at. Shown on the 
same panel, is a stroboscopic surface of section arising 
from the same Hamiltonian, where the green points rep¬ 
resent a chaotic sea while quasi-periodic trajectories are 
shown with gray curves. Note that as expected, the size 
of the chaotic region approximately conforms to the max¬ 
imal phase-space area occupied by the separatrix. In¬ 
deed, this Figure is quite similar to the Poincare surface 
of section shown in panel D of Figure (4). Moreover, the 
differences in the locations of families of quasi-periodic 
trajectories embedded in the chaotic sea, (shown in pan¬ 
els B, C and D of Figure 4), can now be understood as 
changes in the extent of separatrix sweeping that result 
from alteration of the modulation at different values of 

n. 

It is noteworthy that the picture delineated in panel B 
of Figure (5) is not exactly one of adiabatic chaos (Wis¬ 
dom 1985; Henrard & Caranicolas 1990). This is made 
evident by the fact that quasi-periodic resonant trajec¬ 
tories that are not swept by the separatrix do not exist 
in this panel. This is likely due to the fact that the mod¬ 
ulation frequency, a, is so high, that the appearance and 
overlap of secondary resonances acts to wipe out this 
family of orbits (Lichtenberg & Lieberman 1983). For 
reference, panel C of Figure (5) shows a similar strobo¬ 
scopic surface of section where the modulation frequency 
has been artihcially reduced by a factor of 10. 

3.3. Decoherence and Dijfusion 

With a simplified picture at hand, we may now analyti¬ 
cally derive the stochastic properties of the system. First 
and foremost, the above analysis allows us to obtain an 
estimate of the Lyapunov time. Crudely speaking, the 
Lyapunov time can be understood as a characteristic de¬ 
coherence time of the system (Lichtenberg & Lieberman 


1983; Murray & Holman 1997). Accordingly, within the 
framework of Hamiltonian (28), it can be approximated 
as the time interval between successive encounters with 
the separatix, or half the modulation period: 


This simple functional form is in agreement with the es¬ 
timates obtained by Holman & Murray (1996) for the 
Asteroid belt. 

Given that the modulation period is slightly shorter 
than 15 years, the above relation suggests that the Lya¬ 
punov time associated with GJ876 “e” should be of or¬ 
der tl ^ 7 years. To check this estimate numerically, 
we evaluated the Lyapunov time by integrating the full 
linearized variational equations (see Holman & Murray 
1996, Ch. 5 of Morbidelli 2002, Deck et al. 2013) in 
parallel with a direct Wbody simulation®. The numer¬ 
ical calculation yielded a Lyapunov time of tl = 7.26 
years, in excellent agreement with equation (29). As al¬ 
ready mentioned above, this strongly suggests that in 
some similarity with the Kepler-36 system (Deck et al. 
2012), the GJ876 system exhibits rapid dynamical chaos. 

Over timescales longer than a Lyapunov time, it is not 
sensible to treat any one trajectory as being represen¬ 
tative. Instead, a statistical description of irregular tra¬ 
jectories is more appropriate (Lichtenberg & Lieberman 
1983; Murray & Holman 1997). In a uniformly chaotic 
region, the transport of actions obeys the Fokker-Plank 
equation (Wang & Uhlenbeck 1945), which reduces to 
the diffusion equation for Hamiltonian systems (Landau 
1937). Accordingly, the chaotic diffusion coefficient, V, 
quantifies the essential attributes of dynamical evolution. 

In the quasi-linear approximation (Murray et al. 1985), 
the diffusion coefficient is given by the square of the typ¬ 
ical change in action, AZ, that takes place over a deco¬ 
herence (i.e. Lyapunov) time, divided by the Lyapunov 
time. In the non-adiabatic regime (which is evidently 
representative of GJ876), correlations can be neglected 
(see Bruhwiler & Cary 1989 for a discussion) and the 
trajectory can be envisioned to explore the chaotic layer 
uniformly. Accordingly, the average jump in action is of 
order half the size of the chaotic region. We have al¬ 
ready mentioned that within the context of Hamiltonian 
(28), the extent of the chaotic layer approximately cor¬ 
responds to the maximal phase-space area attained by 
the separatrix over a modulation cycle. Thus, a rough 
estimate for AZ is given by: 

(30) 

where Uhyp corresponds to the hyperbolic fixed point of 
the separatrix when at = n (see Figure 5) and is the 
solution to the equilibrium equation: 

0 =/3-f 6[n]u-12[/i]'u ( 77 -/x)-h . (31) 

Combining equations (29), (30) and (31), the explicit 

® The calculation was carried over 10'^ years and the initial tan¬ 
gent vector was randomly oriented. 
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Fig. 6.— Dissipative evolution of the multi-resonant system. Panel A shows the eccentricity of planet “e” as a function of time obtained 
using dissipative A’-body simulations. The behavior of the various curves is indistinguishable from each-other, as the results are not 
sensitive to the value of K. Panel B shows the eccentricity evolution obtained within the framework of a dissipated perturbative model. 
Clearly, a quasi-steady limit cycle, such as that observed within the context of the A-body solution is one possible outcome observed 
in panel B. However, other evolutionary sequences that lead to an approach towards a fixed point are also viable. Panel C shows the 
dissipatively-facilitated evolution of "H. In this panel, the value of "H corresponding to a transition from the limit-cycle attractor regime to 
the fixed-point attractor regime is shown with a dashed line. 


form of the diffusion coefficient becomes: 

^ i^zf a 
tl IGtt 

X (2(6)2/3[/i](H + [h][A]{fi - r,)) + (6)i/y3[/r]2/3 

+ - 16(W + H[A](M-ry))3))2/3)4 

+ 3[/i]2/3)'‘/^(1296[/i]h)"f (32) 

Quantitatively, equation (32) evaluates to I? ^ 3 x 
10“®([n][A]^)/(27r). This suggests that the diffusive 
progress of the eccentricity of planet “e” over an orbital 
period is Ae ^ ^3 x 10“® ^ 10“^. 

We can compare this result with a numerical estimate 
of the diffusion coefficient, calculated on a Poincare sur¬ 
face of section. Specifically, we utilized the perturba¬ 
tive model (22) to compute the average of the square 
of the change in action divided by the time-span be¬ 
tween successive section intersections, using parameters 
corresponding to panel D of Figure (4). This procedure 
yielded I?num = 6.15 x 10“®([n][A]^)/(27r). Thus, the 
analytical estimate underestimates the numerically com¬ 
puted diffusion coefficient by a factor of ^ 2; an accept¬ 
able (and perhaps expected) error given the crudeness of 
the approximation involved in deriving equation (32). 

4. ASSEMBLY OF A CHAOTIC LAPLACE RESONANCE 

As already discussed in the introduction, assembly of 
the “c-b” resonance has been studied extensively in the 
literature (Lee & Peale 2002; Crida et al. 2008 and the 
references therein). Had the discovery of the additional 
planet “e” implied quasi-periodic behavior of the multi¬ 
resonant system (in some similarity to the case of the 
Galilean satellites), the assembly of the Laplace reso¬ 
nance would have been a straightforward extension of 
previous results (Morbidelli et al. 2007). However, the 
chaotic nature of the orbits raises questions regarding 
the specific nature of the natal disk that is both suffi¬ 
ciently laminar to not preclude smooth migration (see 
e.g. Bitsch & Kley 2011) and is simultaneously turbu¬ 
lent enough to prevent the system from settling to quasi- 
periodic depths of the resonant potential well (Adams et 
al. 2008; Rein & Papaloizou 2009). 


For coherence, let us perform our investigation sequen¬ 
tially and first consider only dissipative effects. Accord¬ 
ingly, envisage the “c-b-e” resonant chain embedded in a 
perfectly laminar protoplanetary nebula. Once locked 
in a mean-motion resonance, gap-opening giant plan¬ 
ets (“b” and “c”) tend to carve out a vast mutual gap, 
greatly reducing the disk-facilitated migration rate (Mas- 
set & Snellgrove 2001; Morbidelli & Crida 2007). The 
same argument however does not apply to planet “e”, 
which likely experiences much stronger coupling with the 
disk. Thus, in line with the approximations invoked in 
the previous section, for the purposes of the perturba¬ 
tive model we can take the dynamics of the massive res¬ 
onant pair to be isolated and periodic, while studying 
the dissipative evolution of the outermost planet in the 
test-particle limit. 

4.1. Dissipative Evolution 

Following Lee & Peale (2002), we shall parameterize 
the effects of planet-disk interaction using the following 
simple relationships: 

de e da la , , 

dt Te dt K Te ’ 

where Te is the characteristic eccentricity damping 
timescale, and K is a constant that parameterizes the 
semi-major axis damping timescale in terms of Be¬ 
cause we are interested in understanding the long-term 
evolution of a chaotic multi-resonant state within the 
disk, we shall adopt the present (observed) state as an ad¬ 
equate initial condition. In line with the approximations 
quoted above, we have performed a series of iV-body sim¬ 
ulations where only the outermost planet in the system 
is affected by the fictitious dissipative forces^ (33). 

Recall that the resonant modulation time invoked in 
the previous section is 27r/cr ~ 15 years, much shorter 
than the typically quoted estimates for Tg. Therefore, in 
accordance with adiabatic theory (Henrard 1982), the fi¬ 
nal outcome of our simulations is insensitive to the exact 
value of Tg (which we safely set to Tg = 10^ years). We 
note further that the adopted adiabatic regime® is consis- 

^ Note that this differs from the analysis of e.g. Lee & Peale 
(2002), who initialized the system in a non-resonant state to study 
the capture process. 

® In this limit, Te can be used to replace the effective unit of time 
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tent with the observed orbital state, since resonant cap¬ 
ture probabilities diminish significantly in non-adiabatic 
systems. 

Panel A of Figure (6) shows the eccentricity evolution 
of planet “e” observed in the dissipated A^-body simula¬ 
tions. Specifically, the calculations suggest that follow¬ 
ing an initial transient period of order ~ 5re, the system 
settles onto a quasi-periodic limit cycle. Moreover, the 
evolutionary sequence appears to be independent of the 
assumed value of K, as long as it exceeds K > 5. To 
this end, we note that Lee & Peale (2002) favor a value 
of IK ~ 100 for the natal disk of GJ876 (which we adopt 
for subsequent calculations), rendering our results rather 
robust. 

As with chaotic diffusion, the quoted results of dis¬ 
sipative A^-body simulations can be understood within 
the framework of the perturbative model delineated in 
the previous section. That said, care must be taken in 
implementing equations (33). Explicitly, the definition of 
variables (19) implies that both degrees of freedom will 
be affected by dissipative evolution. However, the pertur¬ 
bative model differs from the A-body system in a crucial 
manner: the test particle approximation employed in the 
former impedes resonant transport of the chain. That is, 
although realistically the application of semi-major axis 
migration to planet “e” affects the radial evolution of the 
whole system, this effect is not captured in the formula¬ 
tion of the Hamiltonian (22)®. Consequently, given that 
the action W is directly proportional to A, a more phys¬ 
ical representation of the real system can be attained by 
only applying dissipative effects to the {u,v) degree of 
freedom. Suitably, the revised equations of motion take 
the form: 


du 

f du^ 

) - - 

dt 

\dt^ 

' n A 

dv 

f dv'' 

) - - 

dt 

\dt, 

n A 


where the subscript T-L signifies the Hamiltonian contribu¬ 
tion. The equations corresponding to the (W, w) degree 
of freedom remain unchanged from (23). 

Using the modified perturbative model, we have sim¬ 
ulated the dissipative evolution of the system, starting 
with different values of H depicted in Figure (4). The 
calculated eccentricity is shown as a function of time in 
panel B of Figure (6), where each curve is labeled by the 
starting value of the Hamiltonian. The corresponding 
evolution of H is shown in panel C of Figure (6). Evi¬ 
dently, two distinct modes of dissipative evolution exist: 
the approach to a limit cycle (as shown by the blue curve, 
T-Li) and the approach towards a fixed point. 

As long as the damping is comparatively slow (as dis¬ 
cussed above), at each point along the evolutionary track, 
a purely Hamiltonian snap-shot of the dynamical por¬ 
trait yields a good representation of the orbital state. 
Thus, the dichotomy inherent to the behavior of the dis¬ 
sipative evolution can be understood by examining the 
surfaces of section depicted in Figure (4). The addition 

associated with dissipative evolution, allowing the obtained results 
to be scaled to other values of t^. 

® An explicit damping of W introduces an unphysical dependence 
on K into the perturbative model. 


of dissipative forces into the model alters the dynamical 
behavior into two ways: it leads to a gradual reduction of 
the phase-space area enclosed by a given orbit in phase- 
space (see e.g. Yoder 1973, 1979; Batygin & Morbidelli 
2011 for a related discussion) and an increase in the value 
of the Hamiltonian. Consequently, if the orbit is initial¬ 
ized somewhere within the chaotic layer (see panels B, C 
and D of Figure 4), as time proceeds the orbit will tend 
to exit the chaotic sea and settle onto the center of the 
corresponding quasi-periodic island. Simultaneously, the 
dynamical portraits will change in such a manner that 
the area enclosed by the constant-energy curves on the 
double section (panel A of Figure 4) will also decrease. 
The two processes are interrelated as the rate of change 
of H grows as the action Z (associated with the {u,v) 
degree of freedom) increases. 

If the starting condition of a "H = T-Li orbit is chosen 
such that it is relatively close to the center of the quasi- 
periodic island depicted in panel B of Figure (4), the 
orbit will rapidly decay onto a limit-cycle that intersects 
the center of the associated Poincare surface of section. 
Concurrently, because such a limit cycle is characterized 
by a consistently low eccentricity, the rate of dissipative 
increase of % will be reduced, rendering such a state of 
the system quasi-steady. An example of such an evolu¬ 
tion is shown by the blue curves in panels B and C of 
Figure (6). The A-body solution shown in panel A of 
Figure (6) also exhibits such behavior. 

If the system is initialized at a higher level of % (e.g. 
7 ^ 3 , 774 , ”^ 5 ), the corresponding Poincare surfaces of sec¬ 
tion depicted in Figure (4) show that the limit cycle 
(which turns out to correspond to a fixed point) resides 
at a high eccentricity meaning that the evolution of 77 
will also proceed at an unhindered rate. Consequently, 
on a timescale not much grater that ^ Tg, the system 
will settle onto an equilibrium described by a maximal 
attainable value of 77 (which also corresponds to a null 
area enclosed by the orbit in the double-section shown 
on Figure 4). Solutions of this kind are shown as pink 
( 773 ), orange (774), and brown {H^) curves on panels B 
and C of Figure (6). 

The value of 77 that separates the two regimes lies be¬ 
tween 772 < 77tr < 773 . Thus, it is possible to have a 
trajectory that first evolves onto a limit cycle, but upon 
crossing the critical value of 77, breaks out of the limit 
cycle^® and begins its approach to the fixed point. An 
example of such an evolution is shown in purple (772) on 
panels B and C of Figure (6). 

Irrespective of the details of the solution, a common 
feature observed in all evolutionary tracks is the ap¬ 
proach to quasi-periodicity on a timescale comparable 
to Tg. Numerical simulations (see e.g. Ogilvie & Lubow 
2003; Papaloizou et al. 2007) suggest that Tg is consid¬ 
erably shorter than the lifetime of a disk. Consequently, 
in absence of additional perturbations, the formation of 
a chaotic Laplace resonance in a perfectly laminar pro¬ 
toplanetary disk appears unlikely. Accordingly, in the 
following discussion we shall invoke turbulent forcing as 
a means of preventing the system from reaching dynam¬ 
ical equilibration. 


This happens because the associated island of stability gets 
engulfed by the chaotic sea. 
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4.2. Damped, Driven Evolution 

Angular momentum transport within protoplanetary 
disks (that facilitates accretion of disk material onto 
the host stars) is typically attributed to turbulence (Ar- 
mitage 2011) that is expected to stem from the magneto- 
rotational instability (Balbus & Hawley 1991) or some 
other process. Turbulent effects can have important con¬ 
sequences for resonant coupling and have generally been 
shown to lead to a breakup of resonant libration (Adams 
et al. 2008; Rein & Papaloizou 2009; Ketchum et al. 
2011). Therefore, provided sufficiently vigorous stochas¬ 
tic diffusion, it is reasonable to expect that turbulent per¬ 
turbations may overcome dissipative interactions, allow¬ 
ing for the formation of a chaotic multi-resonant chain. 

As is with laminar disk-star interactions, the most di¬ 
rect way to calculate the evolution of embedded planets is 
using global 3D magnetohydrodynamic simulations (see 
for example Fromang & Nelson 2006). However, such 
simulations can be computationally expensive and will 
likely preclude a statistically sound exploration of the 
relevant parameter regime. As a result, here we shall 
again opt for a parametrized treatment of stochastic forc¬ 
ing and perform the simulations within the framework of 
the perturbative model. 

The simplest approach to mimicking stochastic forcing 
is by introducing Gaussian white noise, into the equa¬ 
tions of motion (Adams et al. 2008; Rein & Papaloizou 
2009). Following the purely dissipative treatment out¬ 
lined above, we shall only apply these effects to the {u, v) 
degree of freedom. Moreover, as we will seek to examine 
the evolved outcome of these simulations, we shall addi¬ 
tionally introduce a time-dependence, T, to the damping 
and driving terms which will (over a time considerably 
greater than cause these effects to vanish slowly. Ac¬ 
cordingly, the relevant equations of motion become: 


du 

( du'' 


■ 

U 

dt 

\dt, 


- 

^ u 

Te 

dv 

f dv'' 

\ 


V 

dt 

\dtj 


— 

Te 


(35) 


Equations of the form (35) are typically referred to as 
Langevin equations and yield solutions that should be 
interpreted in a statistical sense (Klebaner 2012). It is 
well known that the integration of white noise, yields 
the Weiner process, W (a scaling limit of a random walk). 
With a zero expectation value, the progress of a pure 
Weiner process is ~ VWt, where the ^ is the effective 
diffusion coefficient. ^ is directly related to the variance 
of which we take as an adjustable parameter. 

If the stochastic term in the equations of motion is aug¬ 
mented with dissipative effects (i.e. the terms within the 
brackets of equations 35), the progress of the associated 
quantity becomes limited from above and upon satura¬ 
tion, approaches ~ Consequently, for the system 

at hand, we can define a dimensionless number: 


X = 



(36) 


tions with the nebula^^. Note that the physical meaning 
of X roughly parallels that of the turbulent Schmidt num¬ 
ber: it represents a ratio of turbulent forcing to viscous 
dissipation. Thus, a value of x that leads to the forma¬ 
tion of a system that resembles the real GJ876 resonant 
chain is informative of the properties of the protoplane- 
tay disk in which the system was assembled. 

Of course, the dynamical behavior of the actual 
damped, driven two degree of freedom system is quite 
complex (in part because the value of "H changes sub¬ 
stantially), and it is not clear what value of x leads to 
favorable evolution a-priori. Thus, we have performed a 
series of Monte-Carlo numerical experiments in an effort 
to survey the characteristic outcome as a function of x- 
For all simulations, a starting value of H = Tdi, corre¬ 
sponding to the Poincare surface of section depicted in 
panel B of Figure (4), was adopted and the initial con¬ 
dition was chosen randomly on the section. Indeed, the 
choice of the starting value of ?{ is somewhat arbitrary 
and is generally unimportant because (with the exception 
of a finely tuned set of parameters) the system quickly 
loses memory of its starting state. 

Each integration spanned 30 Tg = 3 x 10"* years and the 
functional form for the time-dependence quoted in equa¬ 
tions (35) was chosen to be T = (exp(—t/f))®, where 
f = 20re. A range of 10“^ < x < I was explored and 
30 realizations of the system were integrated per choice 
of X- Upon completion of the integrations, the end state 
of each simulation was used as an initial condition for 
a purely Hamiltonian integration (equations 23) and a 
surface of section corresponding to the particular orbit 
was constructed. The stochastic properties of the orbit 
were then examined. 

The results of the Monte-Carlo survey are shown in 
Figure (7). Most importantly, panel A shows the prob¬ 
ability of capture into a chaotic Laplace resonance as a 
function of x- It is striking that the probability of repro¬ 
ducing the observed state can be substantial within the 
context of our model, however this remains the case only 
for a rather restricted range of parameters. Particularly, 
the probability of success is nearly 30% for x — 8 x 10“^ 
but drops to zero for x ^ 5 x 10“^ and X ^ 2 x 10“^. 

The reason for a relatively narrow range of x that al¬ 
lows for chaotic resonances as an outcome has to do with 
the evolution of "H. In the case of the over-damped sys¬ 
tem (where x is low), the evolution essentially proceeds 
as described in the preceding discussion of a purely dis¬ 
sipative setup. That is, the trajectory collapses onto a 
nearby quasi-periodic island and 77 evolves to a station¬ 
ary value characterized by an equilibrated orbital state. 
In the case of an over-driven system, a somewhat differ¬ 
ent picture emerges. As turbulent forcing allows the or¬ 
bit to explore the phase-space stochastically, the system 
eventually exits the resonance all together, and settles 
onto an orbit dominated by secular interactions (not ac¬ 
counted for in our model). Such evolutions are marked 
by a substantial decrease in the value of the Hamiltonian. 

The evolution of the eccentricity and the value of 77 
for representative trajectories observed in our simula¬ 
tions are plotted in panels B and C of Figure (7) respec¬ 
tively. The evolutionary track shown in brown, exhibits 


that approximately characterizes the maximal eccentric- Because of this interpretation, it is sensible to limit X by unity 

ity that the orbit would attain exclusively due to interac- above, for bound systems. 
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Fig. 7.— Stochastically driven, dissipative evolution of the multi-resonant chain. Panel A shows the probability of a chaotic end-state as 
a function of the dimensionless parameter x- Panels B and C depict four representative solutions observed in our Monte-Carlo simulations, 
and parallel panels B and C of Figure (6). Specifically, orbits that approach a chaotic end state (blue), a limit cycle (purple), a fixed point 
(green) as well as a non-resonant state (red) are shown. Recall that the calculations are performed in such a way that the stochastic forcing 
as well as dissipative effects diminish with time, leading to a nearly conservative system at the end of the integrations. 


some similarity to that discussed within the context of 
purely dissipative evolution. Specifically, H increases to 
its maximal value and the system settles to the vicinity 
of a fixed point. Note however, that owing to turbulent 
forcing, the system retains a finite libration amplitude 
at the end of the simulation, shown by non-negligible 
oscillation of the eccentricity. The purple evolutionary 
track is also similar to the purely dissipative case, but 
lies in the limit-cycle category. Correspondingly, the ec¬ 
centricity evolution resembles the results of iV-body in¬ 
tegrations, depicted in panel A of Figure (6). 

An over-driven simulation, where the system breaks 
out of resonance and ultimately settles onto a quasi- 
periodic orbit (see also Adams et al. 2008) is shown in 
red. Note that in this simulation, the evolution of the 
Hamiltonian occurs in the opposite sense compared to 
the over-damped simulation. That is, the value of H 
decreases relative to its initial condition. Finally, an 
orbit that settles onto a chaotic Laplace resonance is 
shown in blue. Here, the value of H remains close to 
the Hi — H 3 range where the Poincare surfaces of sec¬ 
tion depicted in Figure (4) show a considerable fraction of 
phase-space that is occupied by chaotic seas. As such, in 
this simulation the system successfully acquires a chaotic 
end-state, facilitated by stochastic elimination of stable 
quasi-periodic trajectories throughout the evolutionary 
sequence. 

5. DISCUSSION 

In this paper, we have examined the dynamical state as 
well as the formation scenario of a multi-resonant chain 
of planets residing in the GJ876 system. In particular, 
we began our investigation by constructing a simplified 
model, based on Hamiltonian perturbation theory, that 
broadly captures the dynamical behavior of the system. 
Upon a detailed examination of this model, we consid¬ 
ered the assembly of the multi-resonant system in pres¬ 
ence of dissipative and stochastic forces. 

With the tally of exoplanetary detections now in ex¬ 
cess of a thousand (Wright et al. 2011; Batalha et al. 
2013), it is tempting to question the value held in de¬ 
tailed examination of the orbital architecture of a par¬ 
ticular system. While such criticism may be appropri¬ 
ate for bodies whose overall state is reminiscent of other 
well-characterized exoplanets, the GJ876 system easily 
escapes such scrutiny as it clearly stands out as a unique 
member of the observational aggregate. To this end, the 


GJ876 “c-b-e” system represents the only known extraso¬ 
lar Laplace-like resonance. Unlike the the Galilean satel¬ 
lites however, the GJ876 multi-resonant system is vigor¬ 
ously chaotic, with a characteristic decoherence time on 
the order of a decade. The primary difference between 
the two configurations stems from the high eccentrici¬ 
ties attained by the massive planets in the GJ876 sys¬ 
tem and the fact that the resonant libration amplitudes 
associated with planet “e” never approached near-null 
values. While interesting in its own right, the nature of 
the GJ876 resonance plays an important additional role 
as a means of placing much-needed constraints on the 
nature of the protoplanetary disk from which the system 
emerged. 

The characterization of the stability, dynamical struc¬ 
ture, and origins of the Laplace resonance inherent to 
the Galilean satellites (within a tangible framework) is 
among the most widely-celebrated achievements of celes¬ 
tial mechanics of the latter half of the 20th century (see 
Peale 1976, 1986 for reviews). As we have shown here, 
the fiercely chaotic Laplace-like resonance of the GJ876 
system can also be understood within the context of a 
simple time-independent perturbative model, character¬ 
ized by two degrees of freedom. Specifically, the mathe¬ 
matical formulation of the governing Hamiltonian is anal¬ 
ogous to a pair of momentum-coupled pendulums. While 
in agreement with direct numerical integration (Rivera 
et al. 2010), this model clearly illuminates the origins of 
stochastic motion and allows one to derive simple ana¬ 
lytic estimates of the Lyapunov time and action diffusion 
coefficient related to the outer-most planet’s eccentric¬ 
ity (Lichtenberg & Lieberman 1983; Murray & Holman 
1997). 

Given the system’s short (~ decadal) timescale for 
stochastic excursions (implied by the analytical estimates 
and observed in the numerical simulations), one may 
readily expect the associated quasi-random signal to be¬ 
come evident in the observational radial velocity time 
series of the system. Specifically, one may naively expect 
that the individual lines observed on the periodogram 
will be chaotically broadened (Laskar 1993b). This con¬ 
sideration may place fundamental limits on any dynam¬ 
ical model’s ability to match the observed signal, even 
provided an outstanding signal to noise ratio (see Laugh- 
lin & Chambers 2001; Rivera et al. 2010; Deck et al. 
2012). In practice, however, the chaotic nature of the or¬ 
bits is probably not the dominant source of error in the 
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data (Baluev 2011) and likely leads to a limited reduction 
in the goodness of fit (Laughlin private communication) . 

A crucial feature that arises naturally within the con¬ 
text of the developed framework is the simultaneous 
rapidity of the chaotic loss of dynamical memory and 
long-term stability. This characteristic contrasts the 
de-populated mean motion resonances (Kirkwood gaps) 
within the Asteroid belt (Wisdom 1983, 1985; Murray 
1986; Henrard & Caranicolas 1990; Morbidelli & Giorgilli 
1990a,b; Nesvorny & Morbidelli 1998) but bears some 
resemblance to the long-term stability of the planets of 
the Solar System (Laskar 1989, 1996; Sussman & Wis¬ 
dom 1988, 1992; Quinn et al. 1991; Murray & Holman 
1999). The analogy is in fact surprisingly robust, es¬ 
pecially for the case of Mercury. That is, much like 
GJ876 “e,” Mercury is characterized by a Lyapunov time 
that is orders of magnitude shorter than the Solar Sys¬ 
tem’s lifetime and has only a slim chance of escaping 
the Solar System within the remaining main-sequence 
lifetime of the Sun (Laskar 2008; Batygin & Laughlin 
2008; Laskar & Gastineau 2009). In a related effort, 
Lithwick & Wu (2011) have shown that Mercury’s sec¬ 
ular evolution can be understood within the context of 
an autonomous Hamiltonian corresponding to a pair of 
momentum-coupled pendulums^^ (see also Sidlichovsky 
1990; Boue et al. 2012), much like the case of the res¬ 
onant evolution of GJ876 “e”. Furthermore, in direct 
correspondence to the stability of GJ876, Batygin et al. 
(2014) have recently shown that Mercury’s long-term sta¬ 
bility also arises from a topological boundary associated 
with the approximate conservation of the Hamiltonian 
itself (see Poincare surfaces of section depicted in Figure 

4). 

The formation of the GJ876 multi-resonant chain al¬ 
most certainly requires convergent approach of the or¬ 
bits facilitated disk-driven migration. While a perfectly 
laminar evolutionary track had previously been invoked 
to explain the (nearly periodic) “c-b” mean motion res¬ 
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Aside from distinct meanings of the variables, a slight dif¬ 
ference between the perturbative representations of Mercury and 
GJ876 “e” is that both pendulums posses the D’Almbert charac¬ 
teristic (Henrard ^ Lemaitre 1983) in the former case, while one 
of the pendulums is simple in the latter case. 


onance (see Lee & Peale 2002 as well numerous other 
references quoted above), our modeling suggests that a 
chaotic configuration is incompatible with formation in a 
purely dissipative nebular environment. Instead, limited 
stochastic perturbations arising from turbulent forcing 
(Adams et al. 2008; Rein & Papaloizou 2009) appear to 
be needed to explain the observed configuration. Ac¬ 
cordingly, we utilized our perturbative model to survey 
a range of assumed disk parameters with an eye towards 
identifying a regime that leads to a successful formation 
of a chaotic multi-resonant chain. This allowed us to sta¬ 
tistically infer (albeit within the framework of an approx¬ 
imate model) the likely combination of relative strengths 
of turbulent perturbations and dissipative effects inher¬ 
ent to the GJ876 protoplanet ary disk. 

The corresponding dimensionless quantity derived in 
section (4.2) is related to the turbulent Schmidt number 
of the disk. Although parameters of this sort are of great 
importance for the quantification of angular momentum 
transport within disks (Armitage 2011), their observa¬ 
tional characterization is at present scarce (Hughes et al. 
2011). Thus, theoretical considerations such as that un¬ 
dertaken in this work are required to inform the relevant 
characteristics. Systems such as GJ876 provide a rare 
opportunity to perform such an analysis. Accordingly, a 
viable and consequential extension of the presented work 
would involve direct magnetohydrodynamic modeling of 
the assembly of the GJ876 multi resonant system within 
its natal disk (see e.g. Nelson 2005; Nelson & Gressel 
2010). Such an investigation would no doubt further 
our understanding of typical formation environments of 
planetary systems and place additional constraints on the 
dominant processes involved in sculpting their orbital ar¬ 
chitectures. 
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